'''
Created on Oct 9, 2009

@author: mkiyer
'''
import os
import sys

from veggie.app.application import CommandLineApplication
import logging

class Novoalign(CommandLineApplication):
    '''
    Novocraft specialises in the development of fast and accurate 
    tools for comparative genomics. Our primary product is an aligner 
    for single ended and paired end reads from the Illumina Genome 
    Analyser. Novoalign finds global optimum alignments using full 
    Needleman-Wunsch algorithm with affine gap penalties whilst 
    performing at the same or better speed than aligners that are 
    limited to two mismatches and no insert or deletes. 

    http://www.novocraft.com
    '''
    _executable = "$SW_ROOT/bioinfo/alignment/novocraft/novocraft-2.05.16/bigdaddy"
    _indexes = {'hg18': 'hg18_contaminant'}
    _indexes_path = "/lab/mkiyer/refdb/novocraft_indexes"

    _quals = {'solexa-quals': 'SLXFQ',
              'solexa1.3-quals': 'ILMFQ',
              'phred33-quals': 'STDFQ'}

    @staticmethod
    def get_indexes_path():
        return Novoalign._indexes_path
    @staticmethod
    def get_index(name):
        return Novoalign._indexes.get(name, None)
    @staticmethod
    def get_quals(fmt):
        return Novoalign._quals.get(fmt, None)

    def run_default(self, fastq_files, quals_format="solexa-quals",
                    genome="hg18", options=None, cwd=None, dryrun=False):
        if options:
            self.addOptions(options)
        if cwd:
            self.cwd = cwd
        # ensure index exists
        index_name = self.get_index(genome)
        if index_name is None:
            raise KeyError('indexes for genome %s not available' % genome)            
        # ensure correct number of sequence files has been passed
        if len(fastq_files) > 2:
            raise ValueError("error in fastq_files param: cannot pass"
                             " more than 2 sequence files to novoalign")        
        # check quals format
        if quals_format not in self._quals.values():
            quals_param = self.get_quals(quals_format)
        assert quals_param is not None        
        # alignment options
        self.addOptions({'-d': os.path.join(self._indexes_path, index_name),
                         '-f': ' '.join(fastq_files),
                         '-F': quals_param})
        # run the program
        logging.debug(self.getCommandLine())
        if not dryrun:
            return(self())

    def run(self, fastq_files, quals_format="solexa-quals",
            output_file=None, genome="hg18", options=None, cwd=None):
        if options:
            self.addOptions(options)
        if cwd:
            self.cwd = cwd            

        # ensure index exists
        index_name = self.get_index(genome)
        if index_name is None:
            raise KeyError('indexes for genome %s not available' % genome)            

        # check fastq files
        paired_end_mode = False
        if len(fastq_files) > 2:
            raise ValueError("error is fastq_files param: cannot pass"
                             " more than 2 sequence files to Novoalign")
        
        # check quals format
        quals_param = self.get_quals(quals_format)

        # alignment options
        self.addOptions({'-d': os.path.join(self._indexes_path, index_name),
                         '-f': ' '.join(fastq_files),
                         '-F': quals_param,
                         '-r': 'All',
                         '-o': 'Native'})
        
        # output path
        if not output_file:
            output_file = "%s_novoalign.txt" % os.path.splitext(os.path.basename(fastq_files[0]))[0]
        output_file = os.path.join(self.cwd, output_file)
        # run the program
        logging.debug(self.getCommandLine())
        retcode, resultpaths = self()